Temperature extremes of 2022 reduced carbon uptake by forests in Europe

The year 2022 saw record breaking temperatures in Europe during both summer and fall. Similar to the recent 2018 drought, close to 30% (3.0 million km2) of the European continent was under severe summer drought. In 2022, the drought was located in central and southeastern Europe, contrasting the Northern-centered 2018 drought. We show, using multiple sets of observations, a reduction of net biospheric carbon uptake in summer (56-62 TgC) over the drought area. Specific sites in France even showed a widespread summertime carbon release by forests, additional to wildfires. Partial compensation (32%) for the decreased carbon uptake due to drought was offered by a warm autumn with prolonged biospheric carbon uptake. The severity of this second drought event in 5 years suggests drought-induced reduced carbon uptake to no longer be exceptional, and important to factor into Europe’s developing plans for net-zero greenhouse gas emissions that rely on carbon uptake by forests.


A Selection of regions
Table S1.Area (million square kilometre), Biomass (GtC) and SPEI for the selected areas.Biomass is taken from Avitabile

B Soil moisture and VPD anomalies
Europe experienced lower than normal soil moisture in 2022 from January onwards in the South, and over the Center and East from March onward.This is indicated by the anomalies shown in Figure S1.However, the soil moisture anomalies did not develop as strongly as in 2018 in Summer, as shown in Figure S2.Whereas JJA of 2018 shows soil moisture anomalies outside the 2σ range in the drought-affected area, soil moisture anomalies outside the 2σ range are shown in 2022 in a smaller area (Fig S2).We note that observations from the SMAP L-band satellite 2 show very similar patterns, with an area-averaged correlation of 0.98, 0.97 and 0.92 for the South, East and North regions, respectively.Contrary to soil moisture, vapor pressure deficit, a measure of atmospheric drought, was much stronger in 2022 than in 2018.In JJA 2022, the VPD was highest over all years considered in 45 / 51 / 15% of the Centre / South / East regions.We calculate VPD from the ERA5 reanalysis data 4 .In JJA 2018, VPD was much less strong, being the highest in only 36 / 24% of the regions North and Centre, respectively.VPD anomalies larger than 1σ were ubiquitous over the European continent in JJA 2022, and VPD anomalies larger than 2σ were found over the Centre, South and Eastern regions (Figure S3).S2.ICOS atmospheric sites used in this study, including their latitude, longitude, station class (M for mountain, T for tall-tower, UP urban/polluted stations, C for coastal sites, S for ocean) and full name of the station 5 .Sites indicated with a * are not included in the CO 2 analyses for a lack of STILT footprints.
Beside CO 2 also carbon monoxide, CO, was measured at some stations.CO is a tracer for incomplete burning, e.g.due to wildfire activity.The large wildfires in France, as described in Section 'Fires', are seen as enhanced CO mole fractions (Figure S4).

C.1 Calculation of background influence
Measured anomalies in atmospheric CO 2 mole fractions are comprised of anomalies in fluxes and in background CO 2 levels, which can either be anomalous themselves or result from anomalous transport patterns in the atmosphere.To account for this, we subtract the anomaly in background CO 2 mole fractions from the observed CO 2 mole fraction anomaly.
Background influences were calculated by transporting a climatological background with the Lagrangian particle dispersion model STILT 6 , driven by the integrated forecast system (IFS, following the IFS cycle development; for more information, see https://www.ecmwf.int/en/publications/ifs-documentation,last access: 15 March 2023), as described in van der Woude et al (2022, referred to as vdW2022 from hereon) 7 .The climatological background was taken from the optimised CO 2 mole fractions from the CarbonTracker Europe (CTE2022) contribution to GCP2022 8 .This climatological (2019-2021) background represents gradients in atmospheric CO 2 mole fractions in the north-south, and east-west directions.As we subtract the linear trend in atmospheric growth rate, the absolute magnitude of the background is removed, and mainly latitudinal differences, combined with differences in wind direction, drive the importance of the background.Table S3 shows that the background is of limited importance and is generally smaller than 0.5 ppm.S2.
Table S3.Relative and absolute importance of the background influence on observed atmospheric CO 2 mole fraction anomalies for used atmospheric stations for selected months in 2022.The range indicates the 5-95 percent quantiles of influence, expressed as percentage of the observed anomaly due to background CO 2 anomalies.The Absolute influence indicates the median of the background influence of all stations, expressed as ppm CO 2 .The mean influence indicates the mean influence of the background relative to transported flux anomalies.The data uses only representative hours over all used stations.The stations are listed in Table S2.

D SiB4 biosphere model
We use a spatially downscaled version of the SiB4 biosphere model 9 , as described in vdW2022.The resulting product has a spatial resolution of 0.1 by 0.2 °and hourly output.Additional to the NEE, as used in vdW2022, we also analyse GPP, TER, soil moisture stress and leaf relative humidity (VPD) stress, as well as the driver data for SiB4 (ERA5 reanalysis data 4 ).Note that SiB4 outputs the net ecosystem productivity (NEP = GPP -TER), whereas EC towers measure the net ecosystem exchange (NEE, which includes fires).As the influence of fires is generally small on EC sites, they can be compared directly and in this work we refer to SiB4 fluxes as NEE to remain consistent with SM2020.We note NEE from the atmospheric perspective, i.e. uptake by plants is negative and positive NEE indicates emissions into the atmosphere.

D.1 Model performance at larger scales
SiB4 shows a good drought response for the 2018 drought, which is demonstrated in SM2020 and vdW2022.Both studies show that transported SiB4 NEE fluxes are in good agreement with tower measurements of CO 2 mole fractions.There is a similar good agreement between the NEE anomalies produced by SiB4 and those from the simple inversions performed in this paper (Section E), indicating that the large-scale absolute response of the biosphere is well represented in SiB4, also during droughts.Furthermore, the spatial patterns of GPP anomalies agree well with those of NIRv for 2022 (see Figures 4 and S7) and with that of SIF for 2018 10 .The spatial correlation between GPP anomalies as calculated by SiB4 and those of NIRv is 0.78 (N=41, P=10 −9 ).We calculated the correlation by first calculating GPP and NIRv anomalies for JJA with respect to 2016-2021.We averaged the anomalies temporally and removed gridcells at which SiB4 simulates a GPP anomaly smaller than 0.5 µmol m −2 s −1 .After this, we binned the NIRv anomalies to 50 bins, removing bins with less than 2 data points.We then correlated the mean NIRv of each bin to the mean GPP of each bin.

D.2 SiB driver data and fluxes at selected EC sites.
We further assessed the SiB4 response by comparing its output to EC measurements at available European forest sites (Table S8).In these comparisons SiB4 is driven by 0.5 by 0.5 degree meteorology from ERA5, and not by the locally observed meteorology at the site, which typically reduces model skill.We nevertheless assess the original gridded SiB4 version here, since we use the model as a tool to look at integral carbon cycle impacts across Europe, and not as a local ecosystem simulator.
We summarised anomalies of SiB4 output at the EC sites, to compare them with the measured anomalies (main Table 1).The driver data anomalies are shown in Table S4 and the flux anomalies, as well as the VPD anomaly and SPEI, are shown in Table S5.
A comparison between the SiB4 anomalies and the EC anomalies shows that SiB4 correctly identifies the years and seasons with extreme events observed in the EC measurements of GPP (see Fig. S5), capturing the sign of summer and autumn NEE anomalies.The slope of the linear regression between latent heat flux anomalies and GPP anomalies at those sites are similar in SiB4 and at measurement towers as well (0.08 and 0.065 µmol m −1 s −1 / W m −2 for EC and SiB4, respectively), indicating a 6/25 Although relative GPP reductions in SiB4 are captured well, as described above, absolute GPP reductions due to drought are indeed better represented in 2018 than in 2022 (see Figure S6).We attribute this to a poorer representation of meteorological drivers, with SPEI calculated from measured data at the EC tower being structurally lower than SPEI calculated from SiB4 driver data.Partly due to this, SiB4 cannot replicate the strong decrease in GPP in the period of the strongest drought such as at Fr-Hes (Fig. S6d) .However, this does not mean that SiB4 has a poor drought response in GPP.Fig. S6 shows the SPEI anomaly of model and EC observations regressed against the corresponding GPP-anomalies, with regression coefficients given in Table S6.It shows that SiB4 captures the slope of ∆GPP/∆SPEI well (SiB4: slope 2.0 vs EC: slope 2.5) and thus the drought response across the SPEI gradient.This, together with its high correlation with NIRV anomalies spatially, lends credence to its drought response at the larger scales we target.Note that because of the offset in the fitted ∆GPP, and the lack of "shutdown" at extreme SPEI's, we view the SiB4 GPP integrals as a lower limit to the estimates in the main text.SiB4 calculates the limitation on GPP from stress due to atmospheric water demand (VPD stress), rootzone soil moisture stress (SM stress) and temperature stress (T stress) 9 .To assess the relative importance of each of these stresses, we calculate the  relative amount of days for which each of the three stresses is limiting GPP. Figure 5 shows the result of this analysis for July and August 2018 and 2022, to assess the difference in dominant stress factor in the different years.
We note that the VPD stress denotes only the leaf evaporative demand, without any temperature dependence, even though temperature influences VPD.Additionally, We note that SiB4 does not simulate the temperature stress to be dominant, even in the southern part of Europe.This is mainly due to the calculation of stress in SiB4, as GPP is limited by the most limiting stress.Although temperature stress is present, the water-related stresses (VPD and SM stress) are more important.The OCN and JULES biosphere models do not match the GPP anomalies in the 2022 warm autumn in some regions over Europe, when compared to the observed NIRv (Figure S9).Where NIRv indicates enhanced GPP over Scandinavia and central Europe (covering France, Germany through Ukraine), Jules simulates a reduced GPP over Scandinavia, SiB4 shows a less strong response, and OCN simulates a reduced GPP over the Eastern part of Europe.Notably, both the GPP response of SiB4 and Jules are similar to the observed NIRv in central Europe.

E Estimates of the NEE anomaly from atmospheric measurements
Using the atmospheric CO 2 mole fraction anomalies described in C, we estimate the total biosphere flux anomaly using: 1) a box model and 2) a simple inversion method, which we describe next.

E.1 Box model
We have first calculated the biosphere CO 2 uptake anomalies using a simple box model, assuming perfect mixing in the atmospheric boundary layer.This simple model directly translates observed CO 2 mole fraction anomalies in the atmosphere to biosphere flux anomalies, and can therefore shed light on the magnitude of the biosphere flux anomalies.This is done using Equation 1.
where ∆CO 2 is the area averaged anomaly in atmospheric CO 2 mole fraction observations at the stations [ppm/month], φ is the area averaged CO 2 flux [mol/m 2 /month] and A represents the area [m 2 ] affected.α [ppm mol −1 ] is a conversion factor for moles of carbon to ppm of CO 2 for a given volume of dry air.We use boxes the size of the four regions, with a mixing layer height of 1000m, which is very similar to the mixing layer height used in STILT (see C.1).We assume that all stations that lie within a region are equally representative of the fluxes only within that region, which is a coarse assumption.We average the observed CO 2 mole fraction anomalies per region to use in the box model.As there are no atmospheric CO 2 mole fraction measurements over the East region, we cannot use the box model to estimate the biosphere anomaly over this region.
To quantify the impact over the drought-affected area, we use the simulated anomaly by SiB4 in this region in the box model.

E.2 Simple inversion
We give a second estimate of the European NEE anomaly informed by atmospheric CO 2 mole fraction measurements using a simple inversion, in which the mismatch between between the simulated and observed CO 2 mole fraction anomalies is minimised.We do this by scaling the NEE anomalies as simulated by SiB4 (See D) in each individual region to be optimally consistent with the observed atmospheric mole fraction anomalies.This optimisation is done using a Kalman filter 13 , which minimises Equation 2: where y 0 are the atmospheric CO 2 mole fraction anomalies (see Section C), with a covariance R. We set R to a diagonal matrix with an uncertainty of 1.5 ppm, which is similar to Munassar et al (2023) 14 .x is a vector of four scaling factors per month (one per region) of the SiB4 NEE anomalies (relative to 2019-2021).H denotes the observation operator, translating the state vector x to atmospheric mole fraction anomalies.We calculated H, the linearised matrix representation of H , by transporting the NEE anomalies from SiB4 to all stations using STILT, driven by IFS.x b is the prior estimate of the state vector, which we set to 0 (no anomaly in NEE), meaning that the calculated anomaly is driven only by atmospheric measurements and SiB4 patterns, but not SiB4 magnitude.x b has a covariance P, which we set to 1.0.This large covariance indicates that little weight is given to the prior estimate of 0 anomaly, resulting in a estimate that is mainly constrained by atmospheric observations.To limit the influence of large biases in observed CO 2 due to fossil fuel plumes, as well as to minimise misrepresented atmospheric transport, we remove observations that have a deviation from the first-guess (SiB4) of more than 3σ .Note that doing a full-fledged inversion, such as in Smith et al (2020, referred to as SM2020 from hereon) 10 , requires detailed (non-climatological) background CO 2 fields, anthropogenic and oceanic CO 2 fluxes.

E.3 Comparison to other inverse estimates of the drought of 2018
We validate our quick inversion with a more complex inversion, as described by SM2020, who use a global atmospheric transport model, as well as influence by fluxes outside our study domain.Moreover, they included anthropogenic and oceanic fluxes and used a lagged ensemble Kalman filter 13 .Their transport was done by TM5 15 at 1x1 degree resolution, whereas we use STILT at 1/12 latitude by 1/8 degree longitude.
Notwithstanding these differences in methodology, we find a very similar response of the European NEE over the droughtaffected area in 2018 (see Section A).This is shown in Figure S10.The figure also shows the anomaly as calculated by SiB4, which shows a remarkable resemblance to both inverse results.Note that the quick inversion starts from a 0 anomaly, and the sign and magnitude are both only informed by atmospheric measurements.The fact that SiB4 and the quick inversion show such resemblance in JJA improves our faith in the ability of SiB4 to simulate NEE anomalies.

F Eddy Covariance data
Eddy-covariance data for 14 forest sites (Table S8) in Europe was used.Most sites experienced decreased precipitation, enhanced incoming radiation and higher air temperatures in summer and autumn 2022, compared to 2016-2021 (excluding 2018).This is shown in Table S9.For all analysis, only day-time data (incoming solar radiation > 10 Wm −2 s −1 ) is used.
Due to the high temperatures in autumn 2022, most of the EC sites experienced higher GPP, but also respiration, compared to the climatology.This is shown in Figure S11.
To calculate SPEI from the EC data, the R SPEI package was used (https://cran.r-project.org/web/packages/SPEI/SPEI.pdf,with potential evapotranspiration calculated using the Hargreaves method.For consistency, the SPEI for the SiB4 simulations at the EC sites was done following the same processing pipeline.For all analysis, we use 3-month SPEI.

G Near-infrared reflectance of vegetation
The near infrared reflectance of vegetation (NIRv) is shown to be a good proxy for leaf phenology, and thereby for GPP 16 .NIRv can give information on the seasonal variation of GPP across a range of ecosystems 17,18 , but more importantly also on drought impacts 10,19 .Because of this relation, we use NIRv to assess spatio-temporal variability of GPP anomalies and the progression of the drought.Hence, we can also use NIRv as an independent validation of the SiB4 GPP response to the drought.

G.1 NIRv-GPP anomalies
From the reduction in observed NIRv, we calculate reductions in GPP over the drought-affected area.For this, we use the linear relation between NIRv and GPP, similar to SM2020.For completeness, the methodology is repeated briefly here.We first calculate the linear regression slope between NIRv and GPP at EC sites for all available site-months.We average these slopes per plant functional type (PFT), to obtain the PFT-specific NIRv-GPP relationship, allowing for an offset (GPP = slope * NIRv + offset) by not forcing this relation through the origin.From the NIRv anomalies, we then calculate the corresponding GPP anomalies, which we integrate over the four regions (see Section A).
To assess the uncertainty on this integrated estimate of GPP anomaly, we calculate the GPP anomaly using four different methods: 1) we use the slopes that we used in SM2020, 2) We used the slopes from SM2020, and linearly detrend the NIRv signal pixel-wise; 3) We used slopes calculated with EC data up to 2023 (see Table S8); 4) We used the slopes from 3) and we linearly detrended NIRv.By calculating the GPP anomaly using these four methods, we obtain an average and a standard deviation.We find that all four methods result in GPP anomalies that differ in the order of 10-20% of their magnitude and are consistent with SM2020.
Over the entire growing season (March-Nov), we find a GPP anomaly of -549.0 ± 64.2 TgC in 2022, compared to an anomaly of -369.3 ± 72.7 TgC in 2018.
A more in-depth analysis of NIRv anomalies is discussed in the following Section.

G.2 NIRv anomalies
NIRv anomalies in 2022, relative to 2001-2021, show significant reductions in NIRv in the Centre, South and Eastern regions (see Figure S13).The figure shows the enhanced uptake in spring 2018 in the Northern region, as well as the reductions in summertime of that year.We find reductions in NIRv in the Centre region that are smaller in 2022 than in 2018, but both years show reductions outside the 2σ range.Contrary to the enhanced uptake in spring 2018 in the Northern region, NIRv retrievals indicate reduced uptake from April onwards, especially in the Eastern and Southern region.Peak reductions in NIRv are observed in July and August in the Centre, East and Southern region, after which a prolonged growing season, indicated by positive NIRv anomalies, starts in from October 2022.Assessing the anomalies per land-use type in the CORINE database 20 , this positive anomaly is especially prominent for broadleaved forests, C3 grasses and C3 crops, with positive anomalies of ± 2σ for each of these land-use types in November over the drought-affected area (see also Figure S13).Contrary, evergreen forests and C4 crops show only small recovery, with anomalies of roughly 1σ .

H Carbon balance impacts
Based on the NIRv-GPP estimate (Section G), SiB4 (Section D) and the simple inversion (Section E.2), we provide a first estimate of the impact of the drought of 2022 on the European carbon budget in summer (Table S10) and autumn Table S11) in the drought-affected area.

I Fires
To analyse the fire contribution to the total CO 2 flux anomaly over Europe, we take fires from the GFAS database 21 , processed following vdW2022.A first comparison between GFAS and GFED 22 shows good agreement between the two, except for over France, where an unresolved issue was found in the GFED data (Guido van der Werf, 2022, personal communication).In contrast, GFAS emissions and emissions estimated by

Figure S1 .
Figure S1.ERA5-Land 3 monthly detrended soil moisture anomalies for January-November 2022 (1-11, respectively), relative to 2000-2021 over the first meter of soil.The contours indicate standard deviations from -2σ to 2σ .The month number is shown in the subtitle of the plot.

Figure S4. 1 -
Figure S4.1-day averaged CO mole fractions measured at four sites (OHP: Obervatoire de Haute Provence, France; CRA: Centre de Recherches Atmosphériques, France; BIS: Biscarrose, France; OXK: Ochsenkopf, Germany) for multiple years (colours).The x-axis shows the month of the year and the day of the month.The sites are listed in TableS2.

7 / 25 Figure
Figure S5.Z-scores (anomaly divided by the standard deviation) for GPP at EC towers and as calculated by SiB4 in different seasons over 2016-2022.

Figure S6 ./ 25 Figure S7 .
Figure S6.a) GPP anomalies in JJA 2022 as function of SPEI for EC (blue) and SiB4 (orange); b) Seasonal cycle of the site FR-Hes.The solid blue (orange) line shows the measured (simulated) GPP in 2022 and the blue (orange) shading the climatology (2016-2021, mean ± 1σ ).We note that the SPEI based on site-measured meteorology is -1.7 and according to SiB4 driver data -0.16 in July of 2022.

Figure S8 .
Figure S8.Total integrated 2022 summertime NEE anomalies in the different regions relative to 2016-2021 as calculated by the biosphere models SiB4, Jules and OCN.Negative numbers signify reduced carbon uptake by the biosphere.

Figure S9 .
Figure S9.Anomalies in October 2022 compared to 2016-2021 for NIRv and GPP as simulated by the biosphere models SiB4, JULES and OCN.Positive anomalies signify higher than average GPP/NIRv.

Figure S10 .
Figure S10.NEE anomalies for the drought-affected area in 2018 (North and Centre) according to the inversion by SM2020, according to SiB4 and according to the quick inversion described here (Simple inversion).The anomalies are relative to 2019-2021 for SiB and Inv, but relative to 2012-2017 for SM2020.

Figure S11 .
Figure S11.GPP, TER and NEE anomalies in autumn (SON) 2022, compared to 2019-2021.The black error bars indicate the standard deviation of the climatology.

Figure S14 .
Figure S14.Ranked NIRv anomalies with respect to 2001-2022 for JJA 2018 and 2022, with lower values indicating lower ranked NIRv observed over that respective gridcell.

Figure S15 .
Figure S15.Same as Figure S14 but for SON.TableS10.Estimated anomalies over JJA 2022 from NIRv (GPP), SiB4 (NEE) and the simple inversion (NEE) in TgC, relative to 2019-2021, and total emissions from wildfires (TgC), according to GFAS.Note that we do not have atmospheric mole fraction observations of CO 2 over the East region, and therefore cannot constrain it in the inversion.Negative GPP anomalies are lower than average values, positive NEE anomalies signify less net carbon uptake than average in 2019-2021.
Vallet et al. (2023) 23 are consistent in France.Anomalously large fires in Europe during the drought of 2022 accounted to ±5.2 TgC emissions, additional to the reduced uptake in the European summer of 2022.CO 2 emissions from fires mainly took place in the South region, where 4.8 TgC was emitted (see also S10.The location of the fires is shown in FigureS16.The figure shows intense fires in the South of France and the Northeast of the Iberian peninsula.Additionally, we use fire counts from the Visible Infrared Imaging Radiometer Suite25 .Cumulative fire counts for Europe are shown in FigureS18.

Figure S16 .
Figure S16.Location and intensity of summer (JJA) 2022 wildfires.The contours indicate the different regions; with the North coloured blue, the centre green, the South yellow and the East coloured red.The size of the marker indicate the maximum intensity, measured as µmol C m −2 s −1

Figure S17 .
Figure S17.Monthly (left) and cumulative (right) total yearly C emissions due to fires over Europe from GFAS 7, 24 .
Figure S19. Figure 5 a and b of the main text in a colorblind friendly colormap.The figures show relative amount of days the stress from VPD, Soil moisture and Temperature is limiting in 2018 (a, b, c) and 2022 (d, e, f), respectively.Lower values indicate more limiting stress.

Figure 5 a
Figure S19. Figure 5 a and b of the main text in a colorblind friendly colormap.The figures show relative amount of days the stress from VPD, Soil moisture and Temperature is limiting in 2018 (a, b, c) and 2022 (d, e, f), respectively.Lower values indicate more limiting stress.

Table S4 .
Similar to S9, but for the high-resolution SiB4, selected at the gridcells the stations are in.

Table S5 .
Similar to Table1, but for the high-resolution SiB4, selected at the gridcells the stations are in.GPP (µmol m −2 s −1 ) NEE (µmol m −2 s −1 ) TER (µmol m −2 s −1 ) SPEI (-) good water-use efficiency according to SiB4.For the 2022 drought, SiB4 tends to underestimate the SPEI levels derived at site level, suggesting that its coarser driver data have dampened the intensity of local impacts recorded in the EC-data, with likely impact on GPP as well.

Table S6 .
Linear regression statistics of SPEI vs ∆GPP corresponding to FigureS6c

Table S8 .
Used eddy covariance (EC) sites and their location in latitude (Lat, degrees north), longitude (Lon, degrees east), altitude above ground level (alt), as well as plant functional type (PFT).*EBF: Evergreen Broadleaf forest Short Name Lat (degree north) Lon (degree east) alt (m.a.g.l.) pft